#1D continuous space model for evolution of range edge with and without the evolution of dispersal
# model gives each individual a Gaussian expectation but allow individuals to vary in their gaussian sd
# plan is to compare situations with and without dispersal evolution and the time to breach the gradient

system(paste("R CMD SHLIB /homes/jc227089/evo-dispersal/SRE/Density1D.c"))
dyn.load("/homes/jc227089/evo-dispersal/SRE/Density1D.so")

# Initialises a matrix with n individuals
# individuals have a location (spX, spY) and are 
# initilised with gene for p.disp
init.inds<-function(n, spX, maxD=20){
	X<-runif(n, 0, spX/2) #position
	D<-runif(n, 0, maxD) # individual sd
	H<-rnorm(n, mean=1, sd=0.05) #Habitat match
	cbind(X, D, H)		
}

# Defines the optimal phenotype and how this varies through space
# hab.val defines the initial optimum
# beta defines the gradient of the shift in optimum
# bf defines the area before optimum begins to shift
hab.space<-function(X, beta, bf, hab.val=1){
	op<-ifelse(X<=bf, hab.val, hab.val+beta*(X-bf))
  op[op<0]<-0
  op
}

# hassel commins population growth where b=1 (contest competition)
# generates expected fecundity
# Might be worth experimenting with a few different growth functions in the end...
hass.comm<-function(lambda, K, N){
	exp<-lambda/(1+(lambda-1)*N/K)
}

# estimates population density at each X as a continuous function using a set bandwidth (an individual's zone of influence)
Nx<-function(X, bw=1){
  dm<-outer(X=X, Y=X, FUN="-")
  dm<-dnorm(dm, mean=0, sd=bw)
  Nx<-apply(dm, MARGIN=1, FUN=sum)
  Nx
}

Nx.prim<-function(X){
	.Call("dens", R_X=as.double(X), R_n=length(X))
}

#determines a survival probability based on distance moved
d.cost<-function(dists, maxDD){
  PS<-1-dists/maxDD
  PS[PS<0]<-0
  PS
}

# reproduces individuals based on density, kills parents
# kills offspring based on habitat
# kills dispersing offspring based on dispersal cost
# disperses individuals based on inherited P
repro.disp<-function(popmatrix, spX, K, lambda, beta, mutn, bf, maxD=20, adap=TRUE, dvec=NULL){
	dens<-Nx.prim(popmatrix[,"X"]) # adult density through space
	hab<-hab.space(X=popmatrix[,"X"], beta, bf) # habitat returned to each individual
	if (adap==TRUE) hab<-1-abs(hab-popmatrix[,"H"]) #decline in survival with habitat fit
  hab<-ifelse(hab<0, 0, hab)
	offs<-rpois(length(popmatrix[,1]), hab*hass.comm(lambda, K, dens)) # surviving offspring
	if (length(offs)==0) return(NULL) #catch extinctions
	inds<-1:length(offs) #vector indexes for offspring
	inds<-rep(inds, times=offs)
	popmatrix<-popmatrix[inds,] #replace parents with offspring. Offspring inherit location and P from parents
	if (class(popmatrix)=="numeric") popmatrix<-matrix(popmatrix, ncol=3, dimnames=list(NULL, c("X", "D", "H"))) #cases where only one individual left
	mut.sample<-rbinom(length(inds), 1, mutn) # individuals copping a mutation
	inds<-which(mut.sample==1)
	if (length(inds)>0 & is.null(dvec)){
		mutation<-rnorm(length(inds), popmatrix[inds, "D"], sd=0.5)
		#mutation[which(mutation>maxD)]<-maxD
		mutation[which(mutation<0)]<-0
		popmatrix[inds, "D"]<-mutation
	}
	mut.sample<-rbinom(length(inds), 1, mutn) # individuals copping a mutation
	inds<-which(mut.sample==1)
	if (length(inds)>0){
		mutation<-rnorm(length(inds), popmatrix[inds, "H"], sd=0.1)
		#mutation[which(mutation>1)]<-1
		#mutation[which(mutation<0)]<-0
		popmatrix[inds, "H"]<-mutation
	}
	if (is.null(dvec)) disp<-rnorm(length(popmatrix[,1]), mean=popmatrix[,"X"], sd=popmatrix[,"D"]) #disperse offspring
	if (!is.null(dvec)) disp<-rnorm(nrow(popmatrix), mean=popmatrix[,"X"], sd=sample(dvec, size=nrow(popmatrix), replace=TRUE)) #disperse offspring
  surv<-rbinom(n=length(disp), size=1, prob=d.cost(abs(disp-popmatrix[,"X"]), maxDD=2*maxD))
	disp[disp<0]<--disp[disp<0]
  disp[disp>spX]<-2*spX-disp[disp>spX]
  popmatrix[,"X"]<-disp
	popmatrix[as.logical(surv),]
}


# plots population density and mean trait values through space
plotter.mean<-function(popmatrix, spX, K, beta, bf, filename=NULL, ...){
	xbin<-cut(popmatrix[,"X"], breaks=spX)
  xmean<-tapply(X=popmatrix[,"X"], INDEX=xbin, FUN=mean)
  disp<-tapply(X=popmatrix[,"D"], INDEX=xbin, FUN=mean)
  hab<-tapply(X=popmatrix[,"H"], INDEX=xbin, FUN=mean)
  if (!is.null(filename)) pdf(file=filename, width=5, height=10)
	par(mfrow=c(3, 1))
  plot(xmean, hab, xlab="Location", ylab="Mean trait value", pch=19, col="darkorange", xlim=c(0, spX), ylim=0:1)
  lines(1:spX, hab.space(X=1:spX, beta, bf, ...))
  plot(xmean, disp, xlab="Location", ylab="Mean dispersal value", pch=19, col="darkgreen", xlim=c(0, spX))
  plot(density(popmatrix[,"X"], bw=1, from=0, to=spX), main="", xlab="Location")
	if (!is.null(filename)) dev.off()
}


